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We present a construction of a matrix product state (MPS) that approximates the largest- 
eigenvalue eigenvector of a transfer matrix T, for the purpose of rapidly performing the infi- 
nite system density matrix renormalization group (DMRG) method applied to two-dimensional 
classical lattice models. We use the fact that the largest-eigenvalue eigenvector of T can be 
approximated by a state vector created from the upper or lower half of a finite size cluster. 
Decomposition of the obtained state vector into the MPS gives a way of extending the MPS, 
at the system size increment process in the infinite system DMRG algorithm. As a result, we 
successfully give the physical interpretation of the product wave function renormalization group 
(PWFRG) method, and obtain its appropriate initial condition. 



§1. Introduction 

Calculation of the maximum eigenvalue of a transfer 
matrix is one of the most fundamental task in statis- 
tical mechanics of lattice models. In two dimensional 
(2D) classical systems, the density matrix renormaliza- 
tion group (DMRG) method is a powerful tool for this 
purpose, 1_4 ) where the maximum eigenvalue is variation- 
ally evaluated by implicit use of the matrix product state 
(MPS). 5 - 12 ) 

In actual applications of the DMRG method to classi- 
cal models, most of the computational time is consumed 
for the determination of the largest-eigenvalue eigenvec- 
tor of the renormalized transfer matrix. 2,3 ' The matrix 
product structure of the variational state helpfully tells 
us how to save the computational time. For instance, in 
the finite system DMRG algorithm, 4 ) a good trial vec- 
tor for the diagonallization of the renormalized transfer 
matrix can be constructed from the MPS. 13 ' 14 '' For the 
infinite system DMRG algorithm, however, such prepa- 
ration of the trial vector is rather difficult, since we have 
to extrapolate the MPS — increase the number of ma- 
trices — toward a lager system size. 

The solution of this problem is partially achieved 
by the product wave function renormalization group 
(PWFRG) method, which increases the number of ma- 
trices in the MPS on the basis of the relation between 
block-spin transformations. 15 -* The numerical efficiency 
of the PWFRG method has been demonstrated for clas- 
sical systems 16, 17 ^ as well as quantum ones. 18 - 22 - 1 In spite 
of the efficiency of the PWFRG method, we have two 
problems to be solved from theoretical view point. One is 
in the physical interpretation of the extrapolation of the 
MPS, in particular, when the system size is small. This 
is because the extension process has been justified only in 



the thermodynamic limit. Another is in the preparation 
of initial matrices, which is necessary to start the recur- 
sive computation of the PWFRG method. So far these 
matrices are set empirically, which may not guarantee 
the numerical stability. 

To find the solution of the above problems of the 
PWFRG method, it is instructive to recall the phys- 
ical substance of the Baxter's corner transfer matrix 
(CTM) method 23 - 1 or its variant, the corner transfer 
matrix renormalization group (CTMRG) method. 24,25 ' 
These CTM type methods calculate the partition func- 
tion of the finite size cluster recursively, where the num- 
ber of iterations corresponds to the linear dimension of 
the system. A preferable point on the methods is that 
the initial condition is well-defined as CTMs with appro- 
priate boundary spin configurations. In this paper, we 
import such well-controllable aspects of the CTM-type 
methods into the extension process of the MPS. 

In the next section we consider a way of the exten- 
sion of the finite size clusters, which can be regarded 
as a variant of the CTMRG method. In §3 we analyze 
the MPS structure of the approximate largest-eigenvalue 
eigenvectors generated by the recursive extension process 
obtained in §2. We successfully give the physical inter- 
pretation of of the PWFRG method, and reconstruct 
its numerical algorithm, which includes the appropriate 
initial condition. In the last section, we conclude the 
obtained results. 

§2. Recursive Construction of State Vectors 

As an example of 2D statistical models, let us consider 
the interaction-round-a-face (IRF) model, which con- 
tains the Ising model as its special case. The IRF model 
is defined by the local Boltzmann Weight W(s'a'\sa) 
assigned for each face of the square lattice, where s, 
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s', a and a' represent g-state spin variables on lattice 
points. (See Fig.l.) For simplicity, wc consider 2-state 
Ising spins (q = 2) and assume the symmetry W (s' a' \ sa) 
= W(a's'\as) and W(s'a'\sa) = W{sa\s'a') throughout 
this article. 26 ) Let us consider the finite size cluster of 
width 2N shown in Fig.l. 27 ) The partition function of 
the cluster is given by the spin configuration sum of the 
product of all the Boltzmann weights in the cluster 

where we have omitted spin variables for simplicity. 
Throughout this section we assume the free boundary 
condition. 




Fig. 1. The arrangement of spins of the local Boltzmann weight 
W(s' a'\sa) (in the left) and the finite size cluster of width 2N 
considered in the partition function Z N in Eq.(2.1) (in the right). 
We show the case N = 4, equivalcntly, 2N = 8. 



We divide the system into the upper and lower halves. 
Figure 2(a) shows the lower half of the system, where 
we label the spins on the cut as s 1 ... s N ajq . . . o~\ from 
the left to right. Considering the product of local Boltz- 
mann weights in the lower half, and taking the config- 
uration sum of spins except for s 1 ...s n o~n we 
obtain a 2 27V -dimensional vector l^^) whose elements 
are *jv[ s i . . . s N a N . . . cj. In the same manner we ob- 
tain the 2 2Ar -dimensional vector (^ N \ that corresponds 
to the upper half of the original cluster. Let us call these 
vectors state vectors in the following. By definitions of 
(& N \ and l^jy), the partition function in Eq.(2.1) is ex- 
pressed by the inner product 

Z N = {* N \* N ). (2.2) 

Let us introduce the row-to-row transfer matrix 

T N [s\ ...s' N cr' N - - - Cil-Si ■■■s N cr N . . . cti] (2.3) 

= VF(s' 1 S2|siS 2 )VF(44l s 2 s 3) ' ■■W(a' 2 a[\a 2 a 1 ) . 

For simplicity, we introduce 'dot product notation' with 
which we abbreviate eq.(2.3) as T N — W ■ W ■ ■ ■ W. 

The state vector \*5f N ) can be used as a good varia- 
tional vector for T N , when the system is off critical and 
N is sufficiently larger than the correlation length. 23 -* 
This is because the variational ratio 

„ = {*n\T n \* n } = (* N \T N \y N ) , A) 

9 (9 N \9 N ) Z N [ ■ ' 

well approximates the largest-eigenvalue of T N in the 
sense that, as N increases, / = — (k In g) / 2N where k 
is the Boltzmann constant approaches to the free energy 



per site in the thermodynamic limit. 




(a) (b) 

Fig. 2. Graphical representations: (a) The state vector \^ N ) 
when N = 4. (b) The reduced density matrix p^. Black dots 
denote spin variables that arc summed. 



We occasionally interpret \f& N ) as a 2^ x 2^ matrix 
$f N , whose elements *S? n (sn ...si\a N •■•c 1 ) are equal 
to ^tvI 5 ! ■ • ■ s N a N ■ • ■ a i\- (Wc have put the vertical line 
| between the first matrix index s N . . .s 1 and the second 
one a N ■ ■ ■o~ 1 , respectively, which corresponds to the left 
and the right half- row spins. We reverse the spin order 
as s N . . . s 1 for the later convenience.) This notation 
is useful when we consider the reduced density matrix 
p^r = Tr L I'J'jvX'I'tvI' where Tr L represents the spin con- 
traction for all the s i in the left side. Figure 2(b) shows 
the graphical representation of p^. 2 ' 23 ' 24 -' Using the ma- 
trix notation ^ N we can simply write pfj as the matrix 
product \E'j v \I/ Ar , since the elements of are explicitly 
given by 

Pn(°'n ■■■o'iWn •■• cr i) = ^n( s n ■■■ s iW'n ••• cr i) 

s N ...s 1 

$ jv(% ■■■ s l\ cr N ■■■<T 1 ). (2.5) 

Similarly we obtain the reduced density matrix for the 
left side p\ = Tr R |* w )(* w | = ^ N ^\. It should be 
noted that all the elements of and p\ are real 

and non negative. The partition function in Eq.(2.2) is 
expressed as 

Z N = Trp^ = Tr/4 = Tr^ N * N . (2.6) 

In the following we consider a recursive construction 
of $> N , starting from TV = 1 up to a certain system size. 
Let us first take the configuration sum for 2 spins on the 
local Boltzmann weight 

*i(*iK) = £W(*Wl*i*i) (2.7) 

to obtain a state vector l^) of the smallest width. (See 
Fig. 2(a) and imagine the case TV = 1.) In order to extend 
the system size and construct \^ 2 ), consider the row-to- 
row transfer matrix T 2 =W -W ■ W, whose elements are 
expressed as 

7 2[ s i s 2 cr 2 cr ll s l s 2 cr 2 cr l] ( 2 - 8 ) 

= W(s[s' 2 \s 1 s 2 )W(s' 2 a' 2 \s 2 a 2 )W(a 2 a' 1 \a 2 c7 1 ) . 
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(See Fig. 3(a).) It is convenient to define half-row trans- 
fer matrices for both the left and the right parts of the 
transfer matrix. For Eq.(2.8), we have 



the reduced density matrix 



T 2 R (a , 2 a' 1 \a 2 a 1 ) = W(a' 2 a' 1 \a 2 a 1 ) 
Tt( S ' 2 s'i\s 2 s 1 ) = W(s , 1 s'2\s 1 s 2 ), 



(2.9) 



where the order of spin indices in T 2 is reversed for con- 
venience, so that T 2 = T 2 holds when W is symmetric. 
To obtain |* 2 ), we "join" T 2 = T 2 L ■ W ■ T 2 R to 1^) as 
follows. First we contract a spin in the half-row transfer 
matrices as 

p 2 VXk 2 ) = £T 2 VXk 2 <7i) 



P2VX|020"l) = £*2( s 2 s ll cr 2 cr l)*2( s 2 s ll cr 2 cr l) 



= E U 2^Wi\Q A 2 fe) U 2 (a 2 a, |&) ,(2.14) 

where we assume the decreasing order for the eigenval- 
ues A 2 (£ 2 ), which are non-negative. As in the DMRG 
method, the matrix U 2 plays a role of the block spin 
transformation from the group of 2 spins a 2 a 1 to a new 
block spin variable £ 2 , which takes 4 states. By oper- 
ating Jj\ to P R , we obtain the renormalized matrix P R 
as 



i£(*2ai|s 2 ) = E^^il^i) . (2.10) ^sVs&ks^) = E U 2^Wl\i2) P fWz (J WlW^2) 



pL 

r 2 



For the later use we define a 2 x 2 matrix P 2 
W ■ P R shown in Fig. 3(b), whose elements are given by 

- P 2[ S 'l s 2 cr 2 cr ll s 2 cr 2] (2- 11 ) 

Then we can obtain \^ 2 ) by applying P 2 to 1^!) as shown 
in Fig. 3(c), where the matrix elements of \E" 2 is explicitly 
given by 

^ 2 (s' 2 s[\a' 2 a[) = J2 ^('WM (2.12) 

X W(S 2 C7 2 |S 2 C7 2 )P 2 R (C7 2 C7 1 |(T 2 )* 1 (S 2 |C7 2 ) . 

It should be noted that the variables contained in i S 1 are 
s 2 and g 2 only. In short, we write Eq.(2.12) as |* 2 ) = 
W-P?\%). 



pL 
r 2 



tm nu 



(a) 



(b) 



(c) 



Fig. 3. Graphical representations: (a) Transfer matrix T 2 = W ■ 
W-W = T 2 L -VK-T 2 R inEq.(2.8). (b) P 2 = P 2 L -VK-P 2 R in Eq.(2.11). 
(c) The extension process of the state vector P^ ■ W ■ P 2 R 1*1^) in 
Eq.(2.12). 



(2-15) 

which is graphically represented in Fig. 4. Here we should 
remark that, in strict sense, the r.h.s. of Eq.(2.15) is 
not a simple matrix multiplication between U\ and P 3 R , 
since U 2 works only on a' 2 and <t[. Keeping this point 
in mind, we compactly write Eq.(2.15) as P R = U 2 P R - 
In the following, we shall put "hat" symbol on matrices 
that play a role of the block spin transformations and 
"tilde" on vectors and matrices generated by the block 
spin transformations. Applying P 3 : 
we obtain an extended state vector 



P 3 L -W-P R to |* 2 ) 



* 3 (4C 2 |^2)= E p 3 L (4C 2 M 



(2.16) 



x W{s 3 a 3 \s 3 a 3 )P 3 (<y 3 ^ 2 \a 3 a 2 )^ 2 (s 3 s 2 \a 3 a 2 ) , 

which contains block spin variables as shown in Fig. 5 (a). 
Since the block- spin transformation is exact, the relation 
z 3 = (*3l*3) = (*3l*3) is satisfied. 



^— 0-0 



Fig. 4. The block spin transformation = t/ 2 P R in Eq.(2.15) 



In the same manner we obtain \*£ 3 ) by applying P 3 = 
P 3 ■ W ■ P R to |^ 2 ), where the matrix elements of P R = 
W ■ P R are given by 

^Vs^Xl^) - W(a 3 a' 2 \a 3 a 2 )P 2 K (a 2 a[\a 2 ) • 

(2.13) 

Hereafter we do not refer to the calculations of Tj^, Pjj, 
and for book keeping purpose. 

Since such recursive extensions of \^ N ) cause the ex- 
ponential blow-up of the vector (or matrix) dimensions, 
we have to introduce the density matrix renormalization 
to reduce the dimensions below a certain number m. 1 ' 2 ) 
The value of m is normally chosen to be 10 ~ 1000. (In 
the following explanation we assume m > 8.) Let us 
consider the case N — 2 as an example. We diagonalizc 



Let us continue the system extension process. We di- 
agonalize the reduced density matrix 

PaVs^l^) = E *3(*3C 2 |^2)*3(*3C 2 k3^) 

= E ^3(^3^1^) A 3 (C 3 ) U 3 (<t 3 Q£ 3 ) (2.17) 

in order to obtain the block spin transformation U 3 that 
maps cr 3 £ 2 to the 8-state block spin £ 3 . This time we 
perform the extension from P R to P R operating U 3 from 
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the left and U 2 from the right 

A R Ke 3 k4%)= E U 3(^2\QW(a' 4 a' 3 \a 4 <J 3 ) 

x PaVafcl^) u *(wM (2-18) 

as shown in Fig. 5(b), or in short Pf = &l(W ■ P 3 K )U 2 . 
(Note that the matrix U 2 in the right hand side contains 
a 3 and a 2 , not a 2 and a 1 .) Up to this stage, we keep all 
the block spin states since we have assumed m > 8. 





(a) 



(b) 



tion Z i = recursively up to arbitrary system 

size i = N. In addition to Z N , we can obtain the nearest 
neighbor spin correlation function 

/„ _ \ _ (*jvI*jv*jvI*jv) 



N u N 1 



at the center of the system. 



(2.23) 



A 706 


' V, 






"^0.704 
C/3 
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'"""•^ 
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m=8 ■ ■ 
m=10 ■-■ 
•••,-80 






0.700 


exact a 







0.002 0.004 0.006 0.008 

1/(2N) 



0.01 



Fig. 5. Graphical representations: (a) Extension from |\E , 2 } to 
|§ 3 ) = P 3 |* 2 ) in Eq.(2.16). (b) Extension from P 3 R to Pf = 
U\{W ■ Pf)U 2 in Eq.(2.18). 



Now we can generalize the extension process up to 
arbitrary system size. If we have j^) up to a certain 
size i > 2, we can construct the extended one 

l* i+ i> = P i+ i\*i) = Pt+i -W-Pf +1 \%) (2.19) 

with 

P? +1 ^U}(W -P^U^, (2.20) 

where the block spin transformations V i and U i _ 1 have 
been obtained by the diagonalizations of pf- and pf"_ 1 , 
respectively. The matrix elements of Eq.(2.19) and (2.20) 
are explicitly given by 

*i+i(*<+iGk<'+i6) = 

E Pi+i( s 'i+iQ \s i+ ifii)W(s' i+1 a' i+1 \s i+1 (7 i+1 ) 
s i+i^i a i+i r ii 

x ^+1 Wi+i^i {s i+ i»i Wi+iVi ) (2.21) 

and 

x ^(^+i^ki+i^J^f(^C i -iki'7i-i)^-i(^'7i-ik) ■ 

(2.22) 

The key point in Eqs.(2.20) and (2.22) is that the label 
of U i is larger than that of U i _ 1 by one, which enables 
us to "extrapolate" the state vector in Eqs.(2.19) and 
(2.21). 

As is performed in the DMRG method, we keep at 
most m numbers of relevant states for the block spin vari- 
ables, and neglect the rest of irrelevant ones in Eq.(2.20). 
We can thus calculate the approximate partition func- 



Fig. 6. Nearest neighbor spin correlation function {s N a N ) in 
Eq.(2.23) of the Ising model at criticality. 



Let us observe the numerical efficiency of the explained 
recursive algorithm in the application to the Ising model. 
Figure 6 shows the calculated (s N a N ) at criticality where 
the correlation length is infinite. The calculated data de- 
viates from the scaling line when the system size exceeds 
the artificial cut-off length introduced by the restriction 
of the degrees of freedom down to m. The observed cut- 
off effect is similar to that of the CTMRG method. 24 - 25) 
This is because the system-size extension process defined 
by Eq.(2.19) and (2.20) is quite similar to the extension 
of the corner transfer matrix in the CTMRG method. 

§3. Product Wave Function Renormalization 
Group 

It is worth looking at the structure of P i = P\ ■ W ■ 
from the view point of the MPS. 6-8 ) Tracing back the 
recursive construction in Eqs.(2.10), (2.15), (2.18), and 
(2.20), we can represent f > i R ( a i£i-il (J i r 7i-i) f° r « > 6 as 

= llUiW ■ U1 2 (W ■ Pt2)Ui-^)Ui-2 (3.1) 

= uUiw ■ ul 2 (w ■(■■■(W- P 2 R ) ■ ■ 0^-3)^-2 . 

Introducing the 2 % by 2 i_1 matrix f^ R (cr^ . . .a[\a i . . . a 2 ), 
which is defined by = W ■ P^ = (W-) l - 2 P 2 R , we 
rewrite Eq.(3.1) as the simpler form 

Pf-sUU — Ulu} Pf U 2 U 3 ---U^ 2 , (3.2) 

where U-_ 1 ■ ■ ■ U 3 U 2 represents the successive block spin 
transformations by U 2 (a 2 a[\^ 2 ), U 3 (a 3 £ 2 \£ 3 ). . . , and 

u i-i( a 'i-iti-2\ti-i)- In the same manner U 2 U 3 ■ ■ ■ U^ 2 
represents the successive block spin transformations by 

^(c^kfe); E^Mlk)-! and ^- 2 (ffi-l1i-2lli-l). 

where one should pay attention to the indices of spin 
variables. The structure of is similar to that of the 



Product Wave Function Renormalization Group: construction from the matrix product point of view 



5 



renormalized half-row transfer matrix 

= iiUiw ■ Uj_ 2 (W ■{■■■{W- T 2 R ) • • 0^-2)^-1 



u 



Ullll (W---W-T?) U 2 U 3 ---U l _ 1 



uU--utu± t? u 2 u 3 ---u i _ 1 , 



(3.3) 



which has been used in the DMRG method. (We 
have introduced T R = W ■ = (VK-) l_2 T 2 R .) 

In Eq.(3.3) the successive block spin transformations 



U. 



■ U 3 Ul are the conjugate of U i _ 1 ■ ■ -U 3 U 2 , and are 
performed by the matrices ^2 C^ ! I £2)1 ^(^^l^)- ■ • 1 
and U i _ 1 (<7 i _ 1 £ >i _ 2 \^ i _ 1 ). Comparing Equations (3.2) 
and (3.3), let us consider the relation between the renor- 
malization process explained in the previous section and 
that in the conventional DMRG method. 

Since the structure of T R in Eq.(3.3) is similar to that 
of Pf in Eq.(3.2), we try to handle Pf- and T R in a 
unified manner. Let us introduce a new matrix V 1 , whose 
elements are given by 



(3.4) 



which is independent of a x , and the right hand side is 
the Kronecker delta. The matrix V 1 has a function of 
taking configuration sum for a 1 as Eq.(2.10), when it is 
operated to the half-row transfer matrices. For example, 
we obtain 



•0-2) 



Tf(a' l ...a' 1 \a l ...a 3 a' 2 W{)V 1 (^\a 2 ) 



Tf{a' i ...a\\a i ...o 2 a x ), 



(3.5) 



which we compactly write as P R = T R V 1 . Then we can 
rewrite P R in Eq.(3.2) as 

P* = Ul l ...Ulul T R V 1 U 2 U 3 ---U i _ 2 . (3.6) 

Note that V x in Eq.(3.4) represents the free boundary 
condition of the system under consideration. For the 
system with fixed boundary condition a x = 1, the r.h.s 
of Eq.(3.4) should be replaced by S(a 2 \a' 2 )6(a 1 |1). 

For the moment let us consider the ideal case m = 
2 J_1 , where we keep all the states for every block spin 
transformations. In such a case the block spin transfor- 
mation is exact, and the relation 



1,1 , r ; ,r,...r ; 



t 



(3.7) 



U 2 U 3 ---U t . 



= L 



is satisfied for the transformation matrices £^2( f7 2 cr il?2)j 

^3(^2^3)- ••. and ^-1(^-1^-21^-1)- The identity 
represents S(a' 1 \a 1 )S((T 2 \(T 2 ) . . . S((r' i _ 1 |cr i _ 1 ), and it 
satisfies = id ■ I { _ 2 = (id-) 1 ^ 1 where id is the local 
identity id(a'\cr) = S(<j'\<t). Putting the identity I i = 
id ■ after T R in Eq.(3.6) and substituting Eq.(3.7), 



we obtain 

P R = UU ■ ■ ■ Ulul T R (J,) 2 V 1 U 2 U 3 ---U i _ 2 
rt . . . jrtirt tR 



u 2 u 3 ■■■u i _ 1 



(3i 



u. 



ulul h v 1 u 2 u 3 ---u i _ 2 



For convenience, let us define a matrix 
Ef = ul 1 ---Uiul k V 1 U 2 U 3 ---U i _ 2 
= id -(UU 1^ V&Ua-- 



■U, 



which satisfies P R 



U 



i-1 



id-E, 



R 



(ul_ 1 Ef_ 1 U i _ 2 ^j , or equivalcntly 
AV&-iM;-i) = *K'h) 



U, 



(3.9) 
id ■ 



E 



tfi-l(*i-l&-2l&-l)£*l(*i-l&-2ki-l»?i-2) 



x ^-2(^-1^-21^-1) ■ 



(3.10) 



By definition the element Ef((r' i £ ii _ 1 \(T i ri i _ 1 ) is zero if 
a[ 7^ a i . In the same manner we define E\ , that is the 
same as P R by the assumed symmetry of the system. 
Substituting Eqs.(3.3) and (3.9) to Eq.(3.8) we obtain 
the simple relation 

P? = T?Ef (3.11) 

among the matrices P/V^i-i l^-i), ^ R K { £i-iK£-_i), 
and Ef'(a'/^' i _ 1 \(T i ri i _ 1 ). Introducing the notations E i = 
E} ■ Ef and f i = f} 1 ■ W ■ f - R we finally reach the relation 



P = r, E t 



(3.12) 



for the ideal case where there is no cut-off in the block 
spin transformation. 

In realistic numerical calculation we have to neglect 
irrelevant states, therefore Eqs.(3.7)-(3.12) do not hold 
exactly, but still P i is well approximated by T i E i if suf- 
ficiently large numbers of the block spin states are kept. 
Under this situation let us consider the operation of P i 
to the state vector 



I*,) = Pi |*i 



(2.19)' 



T i E i \* i _ 1 )=T i \* i ) 



(3.13) 



as shown in Fig. 7(a), where the whole structure of T i 
is shown in Fig. 7(b). We have introduced a new state 
vector 



l*i> = ^l* 



i-il 



(3.14) 



whose structure is shown in Fig. 7(c). If we use the matrix 
representation for the state vectors l^) and we 
can rewrite the above equation by the matrix product 

*i = (^)*i-i(^f) t (3-15) 

among £^ C^Ci-ikC-i), *i-iC s iC;-ih&-i), and. the 
transpose of Ef-(a' i ^' i _ 1 \a i ^ i _ 1 ). The structure of E i = 
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(a) 



(b) 



(c) 



U, 



i-l 



(3.18) 



The initial condition of this relation is obtained from 
the definition of Ef in Eq.(3.9). Substituting Ef = 
UlhVi t0 ^3 = E fU 2 we obtain A 3 = (U^h^i )& 2 
U 2 , and thus we find the initial condition 



UHid-V^ 



A 2 = V, 



(3.19) 



It might be helpful to rewrite Eq.(3.18) by the element: 

A(<7iti-i\Vi) = 



x ^i-ifa-i6-2l»?i-i) Ui-i^iVi-ih) ■ 



(3.20) 



Figure 8 shows the graphical representation of Eqs.(3.18) 
and (3.20). 



Fig. 7. Graphical representations: (a) P. in Eq.(2.19)' in 

case of i = 4. (b) T =fP-W- T\ R . (c) in Eq.(3.14). 

The long rectangular in (a) and (b) represents T- in Eq.(2.3). 
The triangle with number j represents the block spin transforma- 
tion U- or (if. The rectangular area labeled by "W represents 
Small circles and squares are original spins and block 
spins, respectively, where black ones are summed out. Two cir- 
cles that are connected by the double lines represent the same 
spin. We have used Eq. (3.4), and therefore V t is not explicitly 
drawn in the picture. Stacking (b) and (c) wc obtain (a) when 
Eq.(3.7) is satisfied. 



Ef ■ Ef- is understood from Fig. 7(c) without difficulty. 
As we explain in the following, 1^) is easily created from 

l*i-i>- 

It is not required to possess Ef and Ef 1 explicitly for 
the construction of $j in Eq.(3.15). Let us substitute 
the singular value decomposition of to the r.h.s of 
Eq.(3.15). From the assumed symmetry of the de- 

composition is equivalent to the diagonalization of 
as 



Ui-itti-iUli 



(3.16) 



where is a diagonal matrix, and the square of Q i _ 1 
represents the eigenvalues of the density matrix pf_ 1 = 

Ui-i ( fi i-i) 2 ^/-i- Substituting Eq.(3.16) 



to in Eq.(3.15), we obtain 

= (^ L J7 i _ 1 )n i _ 1 (^_ 1 ) t 
= A i Cl i _ 1 A i , 



(3.17) 



where A i = Ef l U i _ 1 and we have used the assumed sym- 
metry Ef = Ef . From the recursion relation of Ef in 
Eq.(3.10) the new matrix A i satisfies the recursion rela- 
tion 



A; = Ef-U, 



uU^d-Eff.U^) 



R 



u. 



U, 



u. 



£>i-2 ^i-1 




Tli-l "Hi 



Fig. 8. Recursive relations of X. X = \u}_ 1 (id - U i _ 1 
for i > 4 in Eqs.(3.18) and (3.20). 



Using the recursion relation Eq.(3.18) and the approx- 
imation from Eq.(2.19)' to Eq.(3.13), let us reformulate 
the recursive creation of the state vectors proposed in 
the previous section. Suppose that we have the renor- 
malized transfer matrix T i _ 1 , the state vector |^ r i _ 1 ) and 
the matrix A i _ 1 . Wc obtain as follows. 

(a) Perform the diagonalization — > U i _ 1 fl i _ 1 U^_ 1 
in Eq.(3.16) to obtain U i _ 1 and 



U t _, by Eq.(3.18). 



(b) Obtain A i = [u}^ (id ■ A i 

(c) Extend the linear dimension of the renormalized 
transfer matrix by ff = (j}^ (w ■ ff_^j in 

Eq.(3.3) and the same for Tf" . 

(d) Obtain $._= A\. (Eq.(3.17).) 

(e) Multiply T t to Admitting the approximation 
from Eq.(2.19)' to Eq.(3.13), substitute the obtained 
vector T i to \^ff) and proceed to the next iter- 
ation. 

With use of this recursive calculation, we obtain the state 
vectors and the partition functions explained in the pre- 
vious section. (The obtained data is of use for the finite 
size scaling at the criticality. 25 ') 

The recursive calculation starts from the initial condi- 
tion given by Eq.(3.19). The calculation can be stopped 
when the observed physical quantities, such as the free 
energy per site, converge to their values at the large 
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system size limit. Alternatively, the difference between 
pfyZj and PjXi/Z i+1 (or in some cases pf'/Z i and 
(^+2/^1+2) can b e use d for the confirmation of the con- 
vergence, since when i is several times larger than the 
correlation length the recursive procedures (a)-(e) be- 
comes just the repetition of the same calculation. It 
should be noted that even in the large i limit the con- 
dition U i+ i = Ui is not always satisfied, as reported by 
Okunishi et al. 2a *> 

The recursive process (a)-(e) is actually a special case 
of the numerical algorithm of the 'product wave function 
renormalization group (PWFRG) method'. 15 ^ In order to 
clarify this point let us consider the numerical problem of 
obtaining the largest-eigenvalue eigenvector of the renor- 
malized transfer matrix T i during the calculation by infi- 
nite system DMRG method. Normally one of the power, 
the Lanczos, the Davidson, or the Arnoldi methods is 
employed for the diagonalization of T i . The computa- 
tional time required for these diagonalization methods is 
greatly dependent on the choice of the initial vector, and 
therefore it is important to create a good initial (or trial) 
vector for the diagonalization of T i . From the arguments 
in the previous section — see Eqs.(2.19) and (2.20) - 
it is apparent that 1^) = P i in Eq.(2.19)' is a 

candidate of the initial vector. (The choice is not so bad 
if the correlation length of the system is much shorter 
than the system size 2i, as was explained in the para- 
graph that contains Eq.(2.4), and expected to be of use 
for modest system size.) Since P i is further approx- 
imated by T i E i , we can say that E i has a function 
of extending the vector dimension of |'I' i _ 1 ) to create 
|i>,-) = E t in Eq.(3.14), and that f i improves j^) 

to give an approximation for l^). We find that l^) 
is also of use as the initial vector, since the diagonaliza- 
tion methods creates % |$ 4 > ~ I*,), (f t ) 2 |$ 2 ) ~ % I*,), 
(Tj) 3 !^) ~ {T i ) 2 |*j), etc., inside their numerical pro- 
cedure, and therefore to start the diagonalization from 
1^) is almost as efficient as to start from 1^). 

Based on the above discussion about the initial vector, 
we find that the recursion relation Eq.(3.18) can be effi- 
ciently used for the numerical acceleration of the infinite 
system DMRG algorithm. The way is simply to replace 
the step (e) of the recursive calculations (a)-(e) by 

(e') Using |<j>j) as the initial vector for the Lanczos di- 
agonalization (or Davidson, Arnoldi diagonalization 
or the power method) 28 ^ of T i . Substitute the ob- 
tained largest-eigenvalue eigenvector of T i to I'fJ 
and proceed to the next iteration. 

The modified recursive method, which consists of (a), 
(b), (c), (d), and (c'), is the numerical algorithm of the 
original version of the PWFRG method. 15 ) So far the 
initial condition of this method — the matrix element of 
A 2 — was set empirically. What we have clarified here is 
that the initial condition A 2 is given by V 1 in Eq.(3.19), 
which represent the boundary condition of the system. 29 ^ 
One can easily verify that the calculation by (a)-(e') gives 
the same result as the infinite system DMRG method, 
since the DMRG method consists of (a), (c), and the 



diagonalization of T i . 

Even with the use of the initial vector <I> i the diago- 
nalization step (e') occasionally consume large amount 
of computational time. We can improve the situation, 
considering the half-way between (e) and (e'). For ex- 
ample, if we employ the Lanczos diagonalization, (e) or 
(e') can be replaced by 

(e") Using as the initial vector for the Lanczos di- 
agonalization of T i , perform only one or only a few 
Lanczos steps. Substitute the obtained Lanczos vec- 
tor to l^&j) and proceed to the next iteration. 

This is a realistic choice if one already has a computa- 
tional program of the infinite system DMRG method, 
and if only the physical quantities at the large system 
size limit i — > 00 is required. 18 - ) 

If we store all the block spin transformations U 2 , U 3 , 
. . ., U i _i obtained successively at the step(a), we implic- 
itly have the variational MPS 6 ~ 10) 



U2---U,. 



(3.21) 



which approximates |\E' i _ 1 ). The formulation can be ob- 
tained by considering the structure of (^ r i _ 1 |T' i _ 1 |^ r i _ 1 ) 
which well approximates ( l J'j_ 1 |T' i _ 1 |\I' i _ 1 ). The fact 
enables us to use the iterative process of the PWFRG 
method (a)-(e'), combine with the finite system DMRG 
method. For example, the process (a)-(e') can be used for 
the initial set-up of the variational MPS, that are further 
improved by the finite system DMRG sweeps. Equally, 
after finishing the finite system DMRG calculation for 
a certain system size 2A r , one can switch to the process 
(a)-(e') again, and after M iterations to obtain a good 
variational MPS for the next finite system DMRG sweep 
at the size 2(N + M). The step (e') can be replaced 
by (e) or (e"), since the finite system DMRG sweep im- 
proves the extended MPS. Such a switching process is of 
use when precise numerical data calculated by the finite 
system DMRG method are required for finite size scaling 
analyses. 

§4. Conclusions and Discussions 

We have explained the physical background of the 
PWFRG method from the view point of the MPS for- 
malism. In section II, we first consider the calculation of 
the partition function of the IRF model on a finite size 
cluster. We observe the fact that the largest-eigenvalue 
eigenvector of a transfer matrix of width 2N can be ap- 
proximated by the state vector that corresponds to lower 
half (in Fig. 2(a)) of the system of width 2N . By use of 
the fact, we construct an iterative numerical method to 
calculate the partition function up to arbitrary system 
size, which is equally of use as the CTMRG method. 24, 25 ^ 
In Section III, we analyze the the whole structure of the 
renormalized transfer matrices, which appears in the for- 
malism shown in §2. Observing the matrix product con- 
struction of the state vector, we find the way of increas- 
ing the mumbcr of matrices contained in MPS when we 
increment the system size. The extended MPS can be 
efficiently used as the initial vector for the diagonaliza- 
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tion of the rcnormalizcd transfer matrix, which appears 
in the infinite system DMRG method. As a result, we 
obtain the numerical algorithm of the PWFRG method. 
The initial condition for the PWFRG method, which has 
been set empirically, is clearly determined by Eq.(3.19) 
that corresponds to the boundary condition of the sys- 
tem. 

The numerical method shown in §2 has clear physical 
meaning as the partition function calculation of finite 
size cluster of square shape. The same result can be ob- 
tained by the PWFRG method with the step (e) in the 
previous section. The PWFRG method with (e ! ) gives 
exactly the same result as the infinite system DMRG 
method, where the calculation of the largest-eigenvalue 
eigenvector of the renormalized transfer matrix is faster 
than the conventional infinite system DMRG method. 
One can use the PWFRG method with (e") if the con- 
vergence of the calculated physical quantities toward the 
large system size limit is slow with the use of (e) or (e'). 
The use of (e" ) accelerates the convergence of the free en- 
ergy per site, but the calculated data before completing 
the convergence do not have clear physical meaning. Af- 
ter the convergence, the PWFRG method with (e), (e') 
and (e" ) give the same result about the thermodynamic 
limit of the system. 

We comment that the PWFRG method can also be 
applied to ID quantum systems, in addition to 2D clas- 
sical systems. 18 ~ 22 ) The simplest example is a uniform 
spin chain, whose Hamiltonian is written as a sum of po- 
sition independent local terms h. In this case we obtain 
the superblock Hamiltonian H i = H\ + h + Hj 1 using 
the recursive relation Hj+ 1 = tjj(h + Hj 1 )^ similar to 
the extension process of the renormalized transfer matrix 
Ti+i = Uj(W ■ Tf-)U i . Implementing the diagonaliza- 
tion of H i into the recursive procedure in the PWFRG 
method (a)-(e') in §3, one can obtain the ground state 
property of the quantum system in the thermodynamic 
limit. As for the classical systems, use of (a)-(e") is effi- 
cient for the acceleration of the convergence toward the 
large system size limit. 

About the recursive creation of variational state for in- 
finite ID quantum Hamiltonians, the CTMRG method 
can also be employed for this purpose. Consider the 
infinite size chess board system generated by the Trot- 
ter decomposition 30 ' 31 ^ applied to the density matrix 
exp(— (3H) at zero temperature. If we divide the 2D 
lattice into 4 CTMs by cuts to the diagonal directions, 
the CTM at the top and the bottom corresponds to the 
variational state, and those at the left and the right corre- 
spond to the exponential of the corner Hamiltonian. This 
approach is quite similar to the system extension method 
explained in §2. In the limit of the zero imaginary time 
step in the Trotter decomposition, one obtains a DMRG 
formalism for the corner Hamiltonian. The formulation 
is similar to the recently proposed renormalization for- 
malism on the corner Hamiltonian by Okunishi, 32 ) but 
the block spin transformation is obtained from the diag- 
onalization of the density matrix, instead of the diago- 
nalization of the corner Hamiltonian. 
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